function corr_idat = fsb_normalize_scans(idat,disp)

% FSB : Function to normalise 4D volumes by average volume intensity
% Determines the average image intensity of every volume and normalizes
% it to give an even overall intensity distribution and thereby reduce
% noise.
%
% EXAMPLE:
% corr_idat = fsb_normalize_scans(idat,disp)
%
% INPUT:
% idat:     4D image data set
% disp:     display graphical information
%
% OUTPUT:
% idat:     a normalized 4D dataset
%
% CALLED BY:
% FSB.m
%
% NOTES:
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
% $ Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,'Scaling volumes to mean volume intensity...');
[a,b,c,d] = size(idat); % use dimensions of idat

roi = idat;

waitbar(10/100)  ;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate average value of brain over time without a for-next loop
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% a1 = int16(sum(roi,1)/a);
% waitbar(20/100)  ;
% a2 = int16(sum(a1,2)/b);
% waitbar(30/100)  ;
% a3 = int16(sum(a2,3)/c);

a1 = sum(roi,1)/a;
waitbar(20/100)  ;
a2 = sum(a1,2)/b;
waitbar(30/100)  ;
a3 = sum(a2,3)/c;
signal_avg = single(squeeze(a3));

waitbar(40/100)  ;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% calculate the corrected values normalized to 1
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
corr = signal_avg/mean(signal_avg);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate a 4-D matrix containing the correction values (faster than
% for-next loops)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
corr2 = repmat(corr,[1 a b c]);
waitbar(50/100)  ;
corr3 = shiftdim(corr2,1);
waitbar(60/100)  ;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Change the original matrix to single to allow for division
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

idat = single(idat);
waitbar(70/100)  ;
corr_idat = idat./corr3;
waitbar(90/100)  ;
%corr_idat = int16(corr_idat);

close(h) ;

if disp == 1;

    roi2= corr_idat;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Calculate average value of brain over time without a for-next loop
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    %     a11 = int16(sum(roi2,1)/a);
    %     a22 = int16(sum(a11,2)/b);
    %     a33 = int16(sum(a22,3)/c);
    a11 = sum(roi2,1)/a;
    a22 = sum(a11,2)/b;
    a33 = sum(a22,3)/c;
    signal_avg2 = single(squeeze(a33));

    % calculate the corrected values normalized to 1
    corro = signal_avg2/avg(signal_avg2);

    figure;subplot(4,1,1);plot (corro);subplot(4,1,2);plot(corr);
    subplot(4,1,3);plot (signal_avg);subplot (4,1,4);plot (signal_avg2);

end

end
